#!/bin/bash -e
# based on methy_anno_170122.sh, has no vcf in input
function info() {
echo Usage: `basename $0` 'methy_result field_chr,field_pos,field_beta'
exit 2
}

while getopts  ":p:f:" opt; do
	case  $opt  in
		p) out_prefix=$OPTARG;;
		f) suffix=$OPTARG;;
		*) info;;
	esac
done
shift $(($OPTIND - 1))


if [ $# -lt 1 ]; then info; fi

. $var

tab2vcf="tab2vcf.sh -p$out_prefix"

methyl_id=methyl_id
beta_id=beta
bed450=$data_path/intervals/1/450.clean.$genome_assembly.bed.gz
vcf_header=$data_path/methy/vcf_header
beta=$out_prefix.beta.txt
# b10k=$data_path/methy/beta/b10k_0.txt
b3k=$data_path/methy/beta/b3k.txt




cat $1|cut -f$2|sed 1d > $beta
bgzip -c $beta > $beta.gz
tabix -f -b2 -e2 $beta.gz

$tab2vcf $beta 1,2

vcf=$out_prefix.convert.vcf

cat $vcf |sed /^##vcfProcessLog/d > $vcf.tmp && mv $vcf $vcf.tmp1 && mv $vcf.tmp $vcf && rm $vcf.tmp1
bgzip -c $vcf > $vcf.gz
tabix -pvcf $vcf.gz

echo methylation id annotate
bcftools annotate -R$bed450 -a$bed450 -h$vcf_header -o$out_prefix.cg450.vcf -Ov -c CHROM,FROM,TO,$methyl_id $vcf.gz

echo beta annotate
bcftools annotate -a$beta.gz -h$vcf_header -o$out_prefix.beta.vcf -Ov -c CHROM,POS,$beta_id $out_prefix.cg450.vcf

echo extract
$java_run/snpsift extractField $out_prefix.beta.vcf $methyl_id $beta_id|sed 1d > $out_prefix.cg450.beta.txt


echo -e '\tbeta\ncancer\tna' |cat - $out_prefix.cg450.beta.txt > $out_prefix.cg450.beta.pd.txt

echo rf pre
rf_pre1.sh $out_prefix.cg450.beta.txt 10000

beta_stat.py $out_prefix.b450k.sub.test.txt $out_prefix.cg450.beta.txt

echo -e '\tbeta\ncancer\tna' |cat - beta_mean > $out_prefix.beta_mean.pd.txt

# rf_methy_011017.py $b10k $out_prefix.cg450.beta.pd.txt 200 4
rf_methy_011017.py $b3k $out_prefix.cg450.beta.pd.txt 200 4
rf_methy_011017.py $b3k $out_prefix.beta_mean.pd.txt 200 4

rf_methy_011017.py $out_prefix.b450k.sub.test.txt $out_prefix.cg450.beta.pd.txt 200 4

rf_methy_011017.py $out_prefix.b450k.sub.test.can_merge.txt $out_prefix.cg450.beta.pd.txt 200 4

rf_methy_011017.py $out_prefix.b450k.sub.test.can_sub.merge.txt $out_prefix.cg450.beta.pd.txt 200 4

rf_methy_011017.py $out_prefix.b450k.sub.test.txt $out_prefix.beta_mean.pd.txt 200 4

rf_methy_011017.py $out_prefix.b450k.sub.test.can_merge.txt $out_prefix.cg450.beta_mean.pd.txt 200 4

rf_methy_011017.py $out_prefix.b450k.sub.test.can_sub.merge.txt $out_prefix.cg450.beta_mean.pd.txt 200 4

bk.sh -tproject/methylation/test/`date +%Y-%m` predict_score.txt

. $cmd_done